Vibrational cooling, heating, and instability in molecular conducting junctions: Full 

counting statistics analysis 
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We study current-induced vibrational cooling, heating, and instability in a donor-acceptor recti- 
fying molecular junction using a full counting statistics approach. In our model, electron-hole pair 
excitations are coupled to a given molecular vibrational mode which is either harmonic or highly 
anharmonic. This mode may be further coupled to a dissipative thermal environment. Adopting 
a master equation approach, we confirm the charge and heat exchange fluctuation theorem in the 
steady-state limit, for both harmonic and anharmonic models. Using simple analytical expressions, 
we calculate the charge current and several measures for the mode effective temperature. At low 
bias, we observe the effect of bias-induced cooling of the vibrational mode. At higher bias, the mode 
effective temperature is higher than the environmental temperature, yet the junction is stable. Be- 
yond that, once the vibrational mode (bias-induced) excitation rate overcomes its relaxation rate, 
instability occurs. We identify regimes of instability as a function of voltage bias and coupling to an 
additional phononic thermal bath. Interestingly, we observe a reentrant behavior where an unstable 
junction can properly behave at a high enough bias. The mechanism for this behavior is discussed. 

PACS numbers: 



I. INTRODUCTION 

Can molecules serve as reliable components in elec- 
tronic circuits? A major obstacle in realizing molecular- 
based electronic devices is junction heating and break- 
down, the result of vibrational excitation by the elec- 
tron current [1-10]. This situation generally occurs once 
the bias voltage exceeds typical molecular vibrational fre- 
quencies and the electronic levels are situated within the 
bias window. If energy dissipation from the conduct- 
ing object to its environment (metals, solvent) is not ef- 
ficient, the molecular conductor experiences significant 
heating, ultimately leading to junction breakdown. A 
related question, the possibility for a nonequilibrium in- 
duced cooling of the junction has been the topic of recent 
experimental and theoretical studies [4, 11-13]. 

In this paper, we study the problem of bias-induced 
molecular cooling, heating, and (potential) junction 
breakdown due to vibrational instabilities, using the 
Donor (D)-Acceptor (A) Aviram-Ratner electronic rec- 
tifier setup [14], see Fig. 1. By coupling electronic tran- 
sitions within the junction to a particular internal molec- 
ular vibrational mode, significant molecular heating can 
take place once the donor level is lifted above the accep- 
tor level, as the excess electronic energy is used to excite 
the vibrational mode. This process may ultimately lead 
to junction instabilities and breakdown [16]. The model 
can also demonstrate current- induced cooling at low bias, 
when tuning the junction's parameters. 

Within this simple system, several issues are of funda- 
mental and practical interest. First, one would like to 
understand the role of mode anharmonicity in the trans- 
port process and in the heating or cooling behavior. Sec- 
ond, the molecular vibration under investigation, the one 
controlling junction stability, can be assumed to be well 



isolated from other modes. Alternatively, this mode may 
be coupled to other phonons, allowing for energy damp- 
ing to a larger environment. These two situations should 
result in distinctive cooling or heating behaviors. These 
issues will be explored here. Other relevant challenges 
which are not considered here are the possibility to se- 
lectively excite vibrational modes in the molecule, using 
voltage bias [17], or more generally, to drive molecular 
motion or trigger chemical dynamics [18]. 
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FIG. 1: Schemes of the two models considered in this work. A 
biased donor-acceptor electronic junction is coupled to either 
(i) a harmonic molecular mode, or (ii) an anharmonic mode, 
represented by a two-state system. In both cases the molecu- 
lar mode may further relax its energy to a phononic thermal 
reservoir, maintained at the temperature T^h- This coupling 
is represented by a dashed arrow. 

Using a full counting statistics (FCS) approach, our 
analysis further contributes to the institution of flue- 



tuation relations in open many-body quantum systems. 
Fluctuation theorems (FT) for entropy production quan- 
tifies the probability of negative entropy generation, mea- 
suring "second law violation" [19, 20]. Such "anomalous" 
processes are relevant at the nanoscale. While originally 
demonstrated in classical systems [21], recent experimen- 
tal efforts are dedicated to explore their validity within 
quantum systems [22]. From the theoretical side, the 
extension of the work and heat FT to the quantum do- 
main has recently attracted significant attention [23, 24]. 
Specifically, a quantum exchange FT, for the transfer of 
charge and energy between two reservoirs maintained at 
different chemical potentials and temperatures, has been 
derived in Ref. [25] using projective measurements, and 
in Ref. [23] based on the unraveling of the quantum 
master equation. It is of interest to test these relations 
in particular cases, e.g., for systems strongly coupled to 
multiple reservoirs, when the reservoirs cooperatively af- 
fect the subsystem [26], including nonmarkovian reser- 
voirs [27-29] , and in models showing coupled charge and 
energy transfer processes, yet the respective fluxes are 
not tightly coupled. The system investigated here corre- 
sponds to the latter case. 

Different flavors of the phonon-assisted-tunneling 
model have been analyzed in the literature [9]. Among 
the various techniques adopted we list solution of the 
dynamics as a scattering problem [30] , extension of the 
basic nonequilibrium Green's function formalism to in- 
clude molecular vibrations [31], or the use of master equa- 
tion approaches [32]. In this paper, we exploit the lat- 
ter method, and present a full-counting statistics of the 
system, allowing for the exploration of charge current, 
energy current and noise processes at the same footing. 
Further, we analytically obtain the cumulant generating 
function (CGF) of the model, allowing for the verification 
of the steady-state charge-energy fluctuation theorem in 
this many-body quantum system. 

The objectives of this work are therefore twofold: (i) 
to analyze a simple model that can elucidate cooling, 
heating and instability mechanisms in molecular recti- 
fiers, specifically, to understand the roles of mode an- 
harmonicity and additional damping routes, and (ii) to 
establish the steady-state entropy production fluctuation 
theorem within a nonequilibrium quantum model, trans- 
ferring charge and energy between the reservoirs in a 
cooperative manner. Recent studies have analyzed the 
role of electron-vibration interaction on the full count- 
ing statistics (FCS) within different approaches [28, 33- 
36]. Complementing these efforts, our treatment offers 
an analytic structure for the CGF, allowing for a clear 
inspection of the microscopic processes involved. 

The plan of the paper is as follows. In Sec. II we 
introduce the D-A molecular rectifier and its two flavors, 
either including a harmonic or an anharmonic internal 
vibration. In Sec. Ill the anharmonic model is analyzed 
within a FCS approach, demonstrating cooling, heating 
and instability dynamics at different parameter regions. 
The case with an additional phonon bath is considered in 



Appendix A. Sec. IV explores the harmonic mode model. 
Sec. V concludes. 



II. MODEL 

Our model includes a biased molecular electronic junc- 
tion and a selected internal vibrational mode which is 
coupled to an electronic transition in the junction. This 
mode possibly interacts with other (reservoir) phonons, 
an extension presented in Appendix A. For a schematic 
representation, see Fig. 1. Generally, this model allows 
one to investigate the exchange of electronic energy with 
molecular (vibrational) heating. The total Hamiltonian 
is given by the following terms. 



H — Hm + Hl + Hr + He + -ff 1, 



H, 



(1) 



The flrst term, Hm, stands for the molecular electronic 
part including two electronic states 



Hm = l^dC^d'^d + <^ac\ca. 



(2) 



Here, cj^, (Q/a) is a fermionic creation (annihilation) 
operator of an electron on the donor or acceptor sites, 
of energies ed,a- The second and third terms in Eq. (1) 
describe the two metals, 7J^, v = L,R^ each including a 
collection of noninteracting electrons 



-ff L = ^ <^ic]^i'^ HR = y^^er 



,t. 



(3) 



IGL 
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The hybridization of the donor state to the left (L) bath, 
and similarly, the coupling of the acceptor site to the 
right (i?) metal, are incorporated into He, 

Hc^'^Vi (cjcd + 4c; j +'^Vr {clca + C^Cr) . (4) 



The Hamiltonian further includes an internal molecular 
vibrational mode of frequency loq. The mode displace- 
ment from equilibrium is coupled to an electron hopping 
in the system with an energy cost k, resulting in heating 
and/or cooling effects. 



Hyib — wofepfeo, 



Hi ^ K c\ca + c\cd (fej + 5o) 



(5) 



Here, h^ (60) represents a bosonic creation (annihilation) 
operator. Note that in our construction the donor and 
acceptor sites are coupled to each other only through 
the interaction with the vibrational mode. We now di- 
agonalizc the electronic part of the Hamiltonian, H^i — 
Hm + Hl + Hr + He, to obtain, separately, the exact 
eigenstates for the L-half and _R-half of Hei , 



He 



E 



eialai 



y^^er-alar- 



(6) 



Assuming that the reservoirs are dense, their new op- 
erators are assigned energies same as those before di- 
agonalization. The donor and acceptor (new) energies 
are assumed to be placed within the band of continuous 
states, excluding the existence of bound states. The old 
operators are related to the exact eigenstates by [37] 



Cd 



I 


[1 


} ^Xrar, 


Cr =} ^rir^r'a,r 



(7) 



where the coefficients, e.g., for the L set, are given by 



A, = 



VLl' 



Su' 



Cd 



VlXl' 



£/ - £;' + 



iS 



(8) 



Similar expressions hold for the R set. It is easy to derive 
the following relation. 



E 



vr, 



ci - ei> + iS 



ppH 



vr, 



(I - Q' 



iTLiei)/2, (9) 



with the hybridization strength ^^{e) — 27r^; vfS{e~ei). 
The expectation values of the exact eigenstates are 

{ajai,) = Si^i'fLiei), (^a^') = 5rylR{<^r), (10) 

where /L(e) — [exp{l3L{€ — fi^)) + 1]~^ denotes the Fermi 
distribution function. An analogous expression holds for 
fni^)- The reservoirs temperatures are denoted by l//3jy ; 
the chemical potentials are ^,^. With the new operators, 
the Hamiltonian (1) can be rewritten as 

Hh = ^ <iia\ai + ^ erolar + Wo^Jfeo 

I r 

+ K^ \\*i\ra\ar + Khalai {bl + bo). (11) 

l,r 

In this form, the model generally describes the process of 
an electron-hole pair excitation by a molecular vibration. 
We denote it by Hh, to highlight the vibrational mode 
harmonicity. A simple version of the model is reached 
by replacing the harmonic mode by a two-state system 
(spin), using the Pauli matrices. 



Ha = y^^eialai 
I 



y erU^Ur 



Wo 



K} A;*Ara[a,. + A*A/a3'.a; 

l.r 



(12) 



The truncated harmonic spectrum imitates an anhar- 
monic mode, as only several (two in the present extreme 
case) states are bounded within the anharmonic poten- 
tial [38]. We denote this Hamiltonian by Ha, to indicate 



on the anharmonicity of the molecular mode. The dy- 
namics of this model should coincide with the behavior 
dictated by Hh, at low temperatures. 

Charge and energy transfer dynamics in these mod- 
els can be followed by studying electronic properties 
[35, 36, 39]. In contrast, here we explore the junction 
response to an applied voltage bias by studying the vi- 
brational mode excitation and relaxation dynamics. The 
analysis of the two-state model, Eq. (12), therefore turns 
out to be simpler than the case when the vibrational 
mode has an infinite spectrum. In what follows, we de- 
rive in details the CGF for the anharmonic-mode case. 
Appendix A generalizes this calculation to include an ad- 
ditional dissipative thermal bath. We then extend these 
results and discuss the model conveyed by Eq. (11). 
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FIG. 2: Scheme of the vibrational mode excitation and re- 
laxation processes. A full circle represents an electron trans- 
ferred; a hollow circle depicts the hole that has been left be- 
hind. 



III. ANHARMONIC-MODE RECTIFIER 
A. Impurity dynamics 

We explore the dynamics of an anharmonic mode, re- 
ferred to as an "impurity", or a two-state-system (TLS), 
within an electronic rectifier, assuming a weak donor- 
acceptor - mode interaction. We rewrite Eq. (12) as 

Ha = -^o-z + a-xFe + ^ eia]ai + ^ eralar, (13) 

I r 

by defining the electron-hole pair generation operator as 
Fe = K^iXiXrajar + A*A;4a;). (14) 

l,r 

The Hamiltonian (13) can be transformed to the spin- 
fermion model [40, 41] of zero energy spacing, using the 
unitary transformation 



C/V,C/ = a„ U^a^U 



CTz, 



(15) 



with U = -j^{<Jx + (Tz). The transformed Hamiltonian 
Ha = U^HaU is given by 

Ha = -^Crx+ CFzFe 



E-^ 



eia\ai + 2^^r<A.o-r- (16) 



In this form, the TLS dynamics can be simulated exactly 
using an influence- functional path integral approach [42] . 
Back to (13), we denote the TLS ground state and ex- 
cited state by |0) and |1), with energies and wq, respec- 
tively. We express the Pauli operators by these states, 
a, = |1)(1| - |0)(0|, a, = |1)(0| + |0)(1|. Next, using the 
quantum Liouville equation, we obtain kinetic-rate equa- 
tions for the states population p„ (n=0,l) [38, 43]. This 
standard derivation involves a second order perturbation 
theory treatment with respect to k, the mode-molecule 
coupling parameter, followed by a Markov approxima- 
tion. The resulting equation for the reduce density ma- 
trix ps take the simple form (fi = 1) 

PS = ~i[v{t),psm 

dTi:TB{[V{t),[V{T),ps{t)pLPR]]} (17) 



Here, V — a^Fe represents the (mode-molecule) coupling 
term in Eq. (13). The operators are written in the in- 
teraction representation, 0{t) = e^iHA-v)tQ^-i(HA-v)t 
and we trace over the electronic degrees of freedom. The 
reservoirs v — L,R are maintained in a grand canoni- 



cal state as p^. 



.(ff,-M.Af.) 



/Z^; Zy is the partition 



function of the v bath. Identifying the diagonal matrix 
elements as population, p„ = \ps\n.m we obtain the ki- 
netic equation 



V\ 



-kt^oPi + ^o^iPo, Pi+Po = 1- 



(18) 



In this model, the off-diagonal elements of the reduced 
density matrix are naturally decoupled from the popula- 
tion dynamics [44]. The excitation {kQ_^i) and relaxation 
(A:f_^o) rate constants are given by Fourier transforms of 
bath correlation functions 



e^"''-(Fe(r)Fe(0))dT 



(19) 



k^ - 

^ — oo 

/oo 
e-^"''^(Fe(r)Fe(0))dr, 
-OO 

enclosing electron-hole pair excitation processes, 
(Fe(i)^e(O)) = «:'TrLTr,,{^^pLPfl 

l.l' r.r' 

X'iXra\{t)ar{t) + \l\ial{t)ai{t) 

A* A,,4(0)a,,(0) + \l,\val,{0)ai,{<d)^ }. 

(20) 



The operators are given in the interaction representation, 
e.g., aj(i) = e*^^*aje~*^^*. As we separately trace over 



the L and i?-baths' degrees of freedom, it can be shown 
that the rate constants can be decomposed into two con- 
tributions. 



satisfying 



k{ 



>0 



== 27rK^ 



X E |AinAr|VL(eO(l - fR{^r))5{0Ja + €i - 6.) 






= 2ttk^ 



Y, \\i?\K\\fL{ei){l - Merm-Luo + ei- e,). 

(22) 



l,r 



Similar relations hold for the right-to-left going excita- 
tions. The energy in the Fermi function /^{e) is mea- 
sured with respect to the (equilibrium) Fermi energy, 
placed at (/x^ -I- pr), and we assume that the bias is ap- 
plied symmetrically, pL — —pR- The four rate constants 
describe distinct electron-hole excitation processes, de- 
picted in Fig. 2. At forward bias, if we set the effective 
density of states (DOS) of the L bath to lie higher in en- 
ergy that the DOS of the right bath, we immediately note 
that the rate k^Z^^ should dominate over fc^J^^f-, poten- 
tially leading to " population inversion" of the vibrational 
mode. Utilizing electronic reservoirs with energy depen- 
dent DOS is thus the basic ingredient of the instability 
formation here, as we show below. For convenience, we 
define the spectral density for the v bath as 



Ue) = 27r«:^|A,f^(e,-e). 



(23) 



jev 



Explicitly, using Eq. (8), we find that this function has 
a Lorentzian lineshape, and that it is centered around 
either the D or A level. 



jR{e) =■ K 



^Lie) 



(e-ed)2+ri(e)2/4 
{e-ear+TRier/A- 



(24) 



Using the spectral density function, we express the terms 
in Eq. (22) as integrals 

1 f°° 

fcf:?o^ = 2^7 .fRi^)[^-fL{e + uJo)]JR{e)JL{uJo + e)de 

1 f'°° 
k^^^ =2^7 .fL{e)[l-.fR{e-uj^)]JL{e)JR{e-uj^)de 

1 f°° 

k^^i = ^J .fR{e)[l-fL{e-iOo)]JR{e)JL{e-u:o)de. 

(25) 



The following 


relations 


hold 


(/3l = 


Mk), 














kt 


¥0 

^R 

¥l 


: e-^^f^e^ 


^0. 

7 


hL^R 
h.R^L 



Pii and A^ = /^l 



— „/3Ap„;3c^o 



(26) 



In equilibrium, detailed bias is therefore maintained. 

The dynamics conveyed by Eqs. (18)-(25) is non- 
separable in terms of the two metals, in contrast to simple 
linear interaction cases [38]. In other words, the reser- 
voirs cooperatively excite or damp energy from the im- 
purity, thus their action is non-additive. 

It should be noted that while we assume a weak inter- 
action limit, between electron- hole pair generation and 
the vibrational mode, our scheme does not enforce weak 
metal-molecule coupling; this part is exactly diagonalized 
to yield the reservoirs spectral function, peaked about 
the D or A levels. If one where to force weak metal- 
molecule interaction, the spectral functions (24) would 
reduce to delta functions, Jl{() = 2'itk5{€ — e^) and 
jR{i) = 27rK(5(e — Co), and the resulting rates would be 
evaluated at the donor and acceptor levels, e.g., fefj^Q^ — 
2iiK'^fL(ed)[l - /i?(ea)]<5(ed - Ca + wo)- This also implies 
that charge and energy currents are not "tightly coupled" 



here, such that for each transferred electron not neces- 
sarily precisely one quanta of energy should be gained 
or drained at either contact. In this aspect, our study 
complements the work reported in [10]. There, using the 
small polaron transformation, the coupling of the molec- 
ular bridge to the leads is assumed to be weak, while its 
coupling to the vibrational mode can be made large. 



B. Resolved charge and energy equations 

We write here a closed expression for the cumulant 
generating function, following the approach developed in 
Refs. [26, 27]. It will allow us to obtain the current, 
its noise power, and to confirm the FTs in this system. 
We define Vtin^ N, w) as the probability that by the time 
t the impurity (TLS) occupies the state n, N electrons 
have been transferred from the L metal to the R side, and 
a net energy uj has been transferred, L to R. Resolving 
Eq. (18) to its charge and energy components, we find 
that this probability satisfies the following equation of 
motion [26, 27], 



Pt{0, N-l,uj-e + wo)/L(e)[l - /^(e - c^o)] JL(e)^K(e - iOo)de 



Pt{0, N+l,iu + e)/fl(e)[l - hie - wq)] Jfl(e) JL(e - iOo)de 
+ / Pt{l,N-l,uj-e-^o).fL(e)[l-.fR{e + LJo)]JL{e)JR(e + LJo)de 



+ I Pt{l,N+l,uj + e)fRie)[l-fL{e + LJo)]JR{e)JLie + LJo)de 

I 



(27) 



One could reason this rate equation as follows. In the first 
equation, the term ■Pt(l,iV, w)fcf^Q stands for the decay 
rate of 7^4(1,^^,0;). The second line describes a process 
where by the time t the TLS occupies the ground state, 
N — 1 excess electrons have arrived at the R terminal, 
and an overall oi to — e + loq energy has been absorbed at 
the R bath. At the time t an electron-hole pair excitation 
generates an electron at the R bath, leaving a hole at the 
L metal. This charge transfer process is accompanied 
by an electronic energy transmission at the amount of 
e — ujo' An electron leaving the L bath has a total energy 
e, however only e — LUQ is gained by the R bath. The rest, 
at the amount of loq, is gained by the vibrational mode. 
A similar reasoning can explain other terms in Eq. (27). 
For convenience, the factor {2tt)^^ in Eq. (25) has been 



absorbed into the definition of J^{uj). 

We Fourier transform the above system of equations 
with respect to both charge and energy, to obtain the 
characteristic function Z{x,ri,t). It depends on the en- 
ergy counting field 77 and the charge counting field x. 



\zix,v,t)) 



00 

N=-c 
00 

N=~c 



aiNx 
aiNx 






(28) 



It satisfies the differential equation 

d\Z{x.r^.t)) 



dt 



W{x,v)\Z{x.V,t)) . (29) 



where the matrix W contains the following elements 





h.L^R I uR^L 




(30) 



Here, 



fH^) 



e"VL(e ± ^o)[l - .fR{^)]JL{e ± ^o) Jfl(e)& 



/OO 
e-^-'il - /L(e ± c^o)]/fl(e) JL(e ± wq) J^le)* 
-CO 

(31) 
The cumulant generating funetion is formally defined as 

G{X,V) = 

lim -In y e'^x / Vt{N,uj)e"^'^dio, (32) 

where we introduced the short notation Vt{N,uj) ~ 
Vt{0,N,uj) +1^1(1, N,uj), that is the probability to trans- 
fer by the time t, N electrons and an energy u from left 
to right, irrespective of the state of the TLS. The charge 
and heat currents can be readily derived, by taking the 
first derivative of the CGF with respect to either 77 or x. 



/j^^^ {N)t ^dG{x,v) 



d{ix) 

IT \ = Ml - d'^ix.v) 



X=0,ri=0 



X=0,r,=0 



(33) 



The quantity (a;)^ denotes the total energy uj transferred 
from L to i? by the (infinitely long) time t; {N)^ similarly 
counts the particles (electrons) transferred in the same 



direction, by that time. The zero frequency noise current 
power density can be similarly obtained. 






{N'),-{Nf, _d^G{x,v) 



t d{zxY 



X=0..)=0 



t 



d{iriY 



X=0,jj=0 



(34) 



The CGF can be expressed in terms of \Z) as 
G(x,??)=lim iln(/|Z(x,r,,t)), 

r— >oo t 



(35) 



with (/| — (11 1, a left vector of unity. It is practically 
given by the negative of the smallest eigenvalue of the 
matrix W, 



Gix.ri) = - 



Wl,l + W2,2 



+ 



\/{wiA -■W2a)'^ + 'iWi^2{X,V)w2s{X,V) 



Wij are the matrix elements of W, see Eq. (30). 



(36) 



C. Fluctuation theorem 

We confirm next the following symmetry 
G(x, v) - G{-x + *(/3lMl - PR^^R), -V + *A/3), (37) 

with A/3 = /3ij — Pl- In order to prove this, we focus on 
the product D{x,'r]) = fi'i^2(x, ?7)w24(X: *?) in Eq- (36), 

Dix,v) = [e^''F-ifj) + e-^^F+iTj)] 

X [e^^F+{rj)+e-^^F2-{n)]. (38) 

Under the transformation x ^^ ^X + iiPbl^L ~ PrI^-r) 
and J] — > —77 + iA/3, using the relation /^(e) = [1 — 
/^(e)]e-'^-(^-^''), we find that 



e'^'F^ir]) -^ e.-'^e-fii^i^i.+fiBMn. / dee-"''e-'^'^'[l-/L(e-a;o)]e-''^('-"''-^^\fe(e)e'^^("~^«)jL(e-wo)Jfl(e) 

J —00 

= e-'^e^'-^'^F^iv)- 

/OO 
dee^^'^e'^'^'^ [1 - /fl(e)] e-'3«('-^«)/L(e + wo)e'^^('^+"«-^^' JL(e + wo)^H(e) 
-00 

= e'^e'^^'^''F+{rj). (39) 



r 



Similarly, one could show that 






(40) 



The extra factors e^^'^^° cancel, and we recover the sym- 
metry 

D{x. V) = D{-x + liht^L - pRm), -V + *A/3), (41) 



confirming Eq. (37). We can now demonstrate the va- 
lidity of a fluctuation relation for this non-equilibrium 
system. The probability to transfer the energy lo by the 
long time i, from L to i?, is given by the inverse Fourier 
transform of Eq. (32), 



P,(^,^)^-^e-^'^ 

— OO 



C(x, r/)e<=('^''')*e-*"''d77, 



(42) 

with limt_5.oo[lnC(x, ?7)]/t — 0. Similarly, the quantity 
Vti—Nj—uj) represents the probability that N charged 
particles and an energy u have been transmitted in the 
opposite direction, right to left, up to time t. Based on 
the symmetry Eq. (37), one can show that [23] 



hm im. ^*(^'-) 



iuAl3 + NipLl^L - /SflMfl) 



t^oo t Vt{-N,-ijj) t 

which is often written in a compact form as 



(43) 



Vt{-N,-u) 



^L0/\[i+N(f3i,ti,L -fial^n) 



(44) 



This expression goes beyond standard metal-molecule 
weak-coupling schemes as the energy and charge trans- 
fer and not tightly coupled, and the energy uj can take 
continuous values, unlike Refs. [39, 46, 47]. 

It should be noted that the above derivation has as- 
sumed charge and energy conservation between the two 
reservoirs. The full particle-energy counting statistics, 
without such an assumption, would begin with the proba- 
bility distribution Vt{n, N^, Nn, ujl, i^r), to find the sys- 
tem at time t in the spin state n = 0, 1, with iV^ electrons 
and Wjy excess energy accumulated at the i^ bath. One 
can readily write an equation of motion for this function, 
analogous to Eq. (27), to be Fourier transformed using 
four counting fields, 

Vt{n, XL, XR, VL,m) = E ^''^'''' E ^''^"''" 



iieTe,'Pt{NL,NR,u}L,u}R) =J2n=o,i'Pt{n,NL,NR,u}L,u}R). 
Enforcing energy and charge conservation, 
N — Nl — ~Nr and cj = lur = —ujl, we recover 
Eq. (44). 



D. Currents, and measures for vibrational cooling, 
heating, or instability 

Currents. Analytical expressions for the charge and en- 
ergy currents are obtained using the definition Eq. (33), 
utilizing Eqs. (30) and (36). These currents are defined 
positive when flowing L to R, and their closed forms are 



(/,) =pi(fcf-J' - fcf-o^) +Po(fco^-^r - fcoVi ), (48) 



and 

(/,) = 

/•OO 

pi I rfwu;/L(w- wo)[l - /fl(w)]JL(cj -u;o)Jfl(w) 

O 

dujuj[l - fL{uj + wo)]/k(i^) Jl(w + ujo)Jji{uj) 

dojujfLioJ + ujq)[1 - fR{uj)]JL{uJ + ujo)Jr{uj) 

dujuj[l - fL{uj - ojQ)]fji{uj)JL{uJ - ujo)Jr{uj) 

(49) 

The TLS population is calculated in the steady-state 
limit. 
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The zero frequency noise current power is given by 
2 



(50) 



(^e 



X / e'^^^-^^duL / e'^^^'^Vt{n,NL,NR,ujL,ujR)dwR. 

(45) 

This quantity satisfles an equation of motion that is anal- 
ogous to Eq. (29). It can be readily proved that the 
negative of the smallest eigenvalue of the corresponding 
matrix W{xL,XR,'riL,r]R) obeys the symmetry 

G{xl,Xr,Vl,Vr) = 

G{-Xl + i/^LfJ-L, -Xr + i/^RfJ-R, -flL + ih, -Vr + ^Pr), 

(46) 

which can be translated into the FT for the probability 
itself, 
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^0^1 + ^1^0 



(h.L^Rh,L^R .l,R^Li,R^L\ ir^\ 



The energy current, directed towards the vibrational 
mode, is zero in the steady-state limit, unless the mode 
is further coupled to a dissipative bath. Formally, it is 
given by the expression 



(/,,b) = -c^oPi [fci^'-^o'' + ^f^'o^] +Po^o [kl;^^ + k^^l 



(52) 



Vt(NL,NR,ujL,ujR) 
rt{-NL,-NR,-LUL,-ujR) 



JNLf!Lf-iL+Nal3Bj^R) 



Measures for vibrational instability. The stability of 
the junction can be estimated, against heating effects, 
by inspecting several measures. First, following Ref. [16] , 
we define the damping rate K^b of the vibrational mode 
as the difference between relaxation and excitation rates. 



X e 



(PrUr + Pi^UJi^) 



(47) 






(53) 



Positive Kyib indicates on the "normal" thermal-like be- 
havior, as relaxation processes overcome excitations. In 
this case, the mode effective temperature (defined below) 
is found to be either below (cooling) or above (heating) 
the environmental temperature, yet the junction remains 
stable in the sense that the ground vibrational state pop- 
ulation is larger than the excited level population. A neg- 
ative value for K^t evinces on the process of an uncon- 
trolled heating of the molecular mode, eventually leading 
to junction instability and breakdown. One can also di- 
rectly inspect the TLS population: population inversion 
reflects on vibrational instability. 

Effective temperature. The TLS population can be fur- 
ther utilized as a measure for the molecular vibration ef- 
fective temperature, l/Peff, (ie/inerf using an equilibrium 
relation. 



Pi 
Pa 



-Pifl^O 



(54) 



A negative value for /3e// attests on population inver- 
sion, thus junction instability. When /3e// is positive, 
one should compare it to the reservoirs' inverse tempera- 
ture/3: If/3e// > /3 the system demonstrates bias-induced 
cooling phenomena. For /3e// < P the vibrational mode 
is heated up relative to its environment. The latter typ- 
ically occurs at an intermediate bias voltage, before in- 
stabilities take place. 
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FIG. 3: (i) Charge current in a rectifying molecular junction. 
Inset: Energies of the donor (full line) and acceptor states 
(dashed line). The dotted lines correspond to the chemical 
potential at the left and right sides, (ii) Damping rate K^u,. 
The junction parameters are r^=0.2, l/Pv = 0.005, k = 0.1, 
LOO = 0.05 and ed(A/i = 0) = -0.2, ea(A/x = 0) = 0.4, all in 
units of [eV]. 



E. Numerical results: isolated mode 

We demonstrate cooling, heating and mode instabil- 
ity upon varying the bias voltage. A generic mechanism 
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FIG. 4: Population of the two-state "vibration" as a function 
of bias voltage. Parameters are the same as in Fig. 3. 



leading to vibrational instabilities (and eventually junc- 
tion rupture) in D-A molecular rectifiers has been dis- 
cussed in Ref. [16]: At large positive bias, when the D 
state is positioned above the acceptor level, electron-hole 
pair excitations by the molecular vibration (TLS here) 
dominate the mode dynamics. This can be schematically 
seen in Fig. 2, where the rate k^Z^^ overcomes other 
rates once the donor spectral function is positioned above 
the acceptor spectral function. As Kyu, becomes nega- 
tive, population inversion is observed. 

The junction setup is displayed in Fig. 3. D and A 
levels are positioned such that in equilibrium, A/i = 0, 
the donor level is placed below the Fermi energy /i, while 
the acceptor level is of a higher energy, ed(A/i = 0) < 
^ < ea{Afi = 0). Under an applied bias, the levels are 
assumed to linearly follow the external potential drive 
(inset) [45]. Therefore, at a particular positive bias the 
levels cross. Beyond that, the levels exchange arrange- 
ment, and the D state is of a higher energy. Throughout 
the paper, the parameters luq, F^, F^^, 1//3, k, ed^a and 
Afi are given in units of eV. 

The junction's current-voltage characteristics is dis- 
played in Fig. 3 (i), manifesting a substantial rectifica- 
tion effect. For negative polarity, Afi = fi^ — (J-r < 0, the 
current is rather small. In contrast, for positive bias the 
current substantially increases once A/x > wq, reaching 
a maximum when the energy levels satisfy ej, — ia ~ wq- 
Level broadening, T^, affects the actual position of the 
maximum. The damping rate, K^b, is displayed in Fig. 
3(ii). It shows the following features: First, for large neg- 
ative bias, A/i < —0.2, Kyu, is negative. This instability 
can be immediately removed, once a very weak coupling 
to a phononic thermal reservoir is turned on, see Figs. 
6 and 11 below. Beyond that, the damping rate Kyib is 
positive between —0.2 < A^ < 0.6, indicating on a sta- 
ble mode of operation. However, for large enough bias, 
A^ > 0.6, once e^ > ea, uncontrolled TLS heating takes 
place, recognized by a sign change in Kyu,. It should be 
noted that the instability takes place in the parameter 
range very relevant to the rectifier operation. It is thus 
important to understand how to tune the system config- 



uration so as to sustain junction functionality. 

Fig. 4 depicts the corresponding population of the two 



levels. At zero bias, /cn 



0, thus the population of 



the excited state is identically zero. At low positive bias 
one finds that kQ_^i < kl_^Q, leading to the "normal" 
situation of po > pi. However, once the bias is large 
and the donor state is positioned above the acceptor site, 
(A/i ^ 0.6) the excitation rate fcg^j^ exceeds the relax- 
ation rate kf^Q and population inversion takes place. We 
note that for a negative bias, small population inversion 
is also observed, as electrons damp energy to the TLS 
when crossing the junction. However, since (/g) is rather 
small (Fig. 3), we do not expect molecular instability in 
this regime, see also Fig. 11. 

The details of the damping rate Kyib depend on the 
level broadening and the reservoirs temperature as we 
show in Fig. 5. The position of the turnover, between 
positive to negative damping, appears at a similar value 
for the bias, and it is generally independent of the reser- 
voirs temperatures and T^. However, the width of the 
curve largely depends on these parameters. 

It should be noted that the development of the insta- 
bility, as reported in Figs. 3, 4 and 5, does not depend 
on the concrete value of k, the strength of the molecule- 
mode coupling, and the behavior persists in the limit of 
vanishing vibronic coupling, k — >■ 0. In the next sec- 
tion we allow the vibrational mode to thermalize with a 
phononic environment at a rate Tph- In this case, the 
competition between k and Tph determines the onset of 
instability, see Eq. (58). 
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FIG. 5: Damping rate in a rectifying junction for differ- 
ent broadening parameters, r^=0.2, fii, = 200 (full) r^=0.4, 
I3„ = 200 (dashed) r^=0.2, /3,^ = 5 (dashed-dotted). Other 
parameters are the same as in Fig. 3. 



F. Numerical results: dissipative mode 

Up to this point, we have assumed that the molec- 
ular vibrational mode (TLS here) is well isolated from 
other vibrations. In reality, internal modes typically ex- 
change energy with "secondary" reservoirs modes, either 
internal, or part of a larger environment, opening up an 



additional route for energy dissipation. It is expected 
that in the presence of such a thermal bath, the region of 
vibrational instability {Kyib < 0) would become limited. 
A simple model that is capable of describing a hier- 
archy of energy transfer processes, electronic energy — > 
specific vibrational excitation — >■ thermal bath, is given 
by an extension of the model (13), 



Ha+i 






F,) 



E' 



eiajai + ^ e^^a^ + ^ coabiba- (55) 



The notation ^^Ha+b^^ indicates that the anharmonic 
mode is coupled to a thermal bath (B). The operator F^ 
describes electron- hole pair excitations as in Eq. (14). 
The thermal bath operator, coupled to the TLS transi- 
tions, includes displacements of reservoir modes. 



F, 






Vaibi + 6q). 



(56) 



with &Jj (6q) as a bosonic creation (annihilation) operator 
for the a phonon-reservoir mode. 

Derivation of the full counting statistics can be reit- 
erated, while including energy dissipation from the TLS 
to the phonon bath. For details, see Appendix A. We 
find that the expression for the charge current stays in- 
tact, satisfying the formal expression (48). However, the 
steady-state populations are corrected by a phonon re- 
laxation rate constant as 



Pi 



^O-fl 



Tphi^o)nphi^o) 



+ kf^o + Tph{uJo)['^nph{uja) + 1] 



(57) 



^O-s-l 



The electronic transition induced rates fc^_j.„/ are those 
defined in Eq. (19); the phononic relaxation rate con- 
stant is Tph{uj) = 27r^^ u^(5(a;Q — w). The function 
nph{uj) = [e'^'''''^ — l]^^ stands for the Bose-Einstein distri- 
bution with /3ph as the temperature of the phonon bath. 
Fig. 6 presents the steady-state population for two 
choices of Tph- When this parameter is small, popu- 
lation inversion still takes place around donor-acceptor 
level crossing. However, the phenomenon disappears at 
large enough bias. Thus, quite interestingly, the domain 
of instability extends intermediate bias values, while the 
system becomes stable again at very high bias. This can 
be reasoned by inspecting Kmt. It is defined as the dif- 
ference between TLS relaxation and excitation rate con- 
stants. In the presence of a thermal bath it is given by 



Kvib = {kl^Q+Tphinph + l)) - {kQ^-f^+Tphnph) . 

fco^i+Tp/^. (58) 



= kl 



For convenience, the uq dependence of the rates is left 
out. While fci_^o < '^o^i ^i^Y ^^old at large bias, both 
these rates diminish with A/i, and the net damping rate 
can become positive due to the Fp^ contribution. For 
large enough Tph, instability does not take place at any 
voltage. 
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FIG. 6: TLS population as a function of bias voltage for V^h ~ 
0.001 (narrow lines) and Tph = 0.1 (heavy lines). The excited 
(ground) state population is presented by dashed (full) lines. 
Other parameters are the same as in Fig. 3 with fiph = 200. 
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absence of a phonon thermal bath we find that at zero 
bias voltage the molecular mode is thermalized at the 
metals' temperature, /3e// = /3. This effective inverse 
temperatm'e quickly drops with increasing bias, becom- 
ing negative around the value of A/i — 0.6, where the 
D-A levels cross. This behavior indicates on the instabil- 
ity of the junction from that point. Next, we weakly cou- 
ple {Tph— 0-001) the single mode to an additional thermal 
bath maintained at fiph = 40. The following observations 
can be made: (i) For negative bias, the mode is close to be 
equilibrated with the phonon bath, as electron-hole exci- 
tations are sparse, (ii) In accordance with Fig. 6, /3e// 
can reach a negative (unstable) value around A/i = 0.6. 
However, /3e// becomes positive at large enough bias, in- 
dicating that the system re-enter a stability region, (iii) 
At low bias, —0.05 < A/i < 0.1, the mode temperature 
is lower than its phononic environment, as /3e// > Pph- 
(iv) At strong mode-thermal bath coupling, Tph = 0.4, 
the mode is closed to be thermalized with /3ph at all bi- 
ases. 

We now demonstrate mode cooling, to a temperature 
below the phonon bath and metals temperature. Keep- 
ing both electron and phonon reservoirs at a fixed tem- 
perature of /3 = 40, the temperature difference AT = 
Teff — Tph is presented in Fig. 8, for various frequencies 
and Tph values. Generally, we note that at low positive 
bias, A/i < 0.1, one may cool the mode by 40 K, the 
result of its coupling to a nonequilibrium environment. 
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FIG. 7: Effective TLS temperature, Fph — (dashed line); 
Tph = 0.001 and /3ph = 40 (full hne) and Tph = 0.4 and 
Pph = 40 (dashed-dotted line). The inset zooms on the latter 
two cases. The dotted lines mark the values Pph = 40 and 
/? = 0. Other junction parameters are the same as in Fig. 3, 
with /3„ = 200. 



IV. HARMONIC-MODE RECTIFIER 

We study next the dynamics of model (11), assuming 
a harmonic mode coupled to the electronic system. The 
relevant equations of motion for the mode levels popula- 
tion are [381 
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FIG. 8: Cooling of the molecular vibration for wq = 0.05 
(dotted line), ujq — 0.15 (dashed line), ujq — 0.3 (full line), 
(a): Tph = 0. (b): Tph = 0.001. Other junction parameters 
are Tl ~ Tr = 0.1, and Pph = Pu = 40. The levels are shifted 
with the bias voltage as depicted in Fig. 3. 

The effective TLS temperature, defined in Eq. (54), 
is displayed in Fig. 7 for several cases. First, in the 



Pn = -[nkl + {n + l)kl]pn 

+ [n+ l)klpn+i + nklpn-i- 



(59) 



Here, the decay rate constant is independent of the level 
index A;^ = k\^Q, and similarly, fc^ = k^^-^, defined 
in Eq. (19). In order to calculate the CGF, we de- 
fine Vt{n,N,u)) as the probability that by the time t 
the harmonic mode occupies level n, N electrons have 
been added to the R bath and an additional energy uj 
has been acquired by the R bath. This quantity fol- 
lows a differential equation analogous to Eq. (27). The 
characteristic function is an array whose nth element is 
\Z{x,V.t))n = Y.N^''''' !^oo^t{n,N,^)e-^^ du. It sat- 
isfies a differential equation corresponding to Eq. (29) 



d\Z(x,V,t)) 
dt 



= -W(x,r/)|Z(x,7],i)), 



(60) 



with the n X n matrix W(x, ?/), 
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(61) 



We can readily confirm the fiuctuation theorem, by in- 
specting the eigenvalues of det [XI — W] . For convenience, 
we define the auxihary matrix A = XI — W. Since it is 
tridiagonal, its determinant can be evaluated in a recur- 
sive manner as 

det[A]i,...,„ = a„,„det[^]i,...,„_i 

- a„^„_ia„_i^„det[^]i,...,„_2 (62) 

where [AJi^..^^; denotes the submatrix constructed by the 
first k rows and columns of A. Thus, the symmetry 
of det [A] with respect to x and rj is determined by the 



symmetry of the products a„.„_ia„_i_„ 
with Wi o the matrix elements of W, 



^n,n— l^n— l,n; 



dniXiV) = Wn,n-l{x,V)Wn-l,n{x,V) « 

[e'>^Ff (ry) + e-'>^F+{7j)] [e'^F+{7j) + e-'^F,-iij)] . 

(63) 

Using the relations (39)-(40), wc conclude that 

dnix, V) == dn(-x + i{hl^L - I^RfiR), -T] + iA/3). (64) 

Given the recursive nature of det [A], this symmetry holds 
for all the eigenvalues of VV, confirming the fluctuation 
theorem (37). We now obtain the steady-state popula- 
tion of the harmonic mode, by solving Eq. (59) in the 
long time limit, p„ = 0. This results in [38] 
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0,1..., oo (65) 



(66) 



if fc^ < fc^. In the opposite limit, the system passes into 
the unstable regime, and the levels' population diverges. 
In that sense, the harmonic model is unphysical as the 
number of states is not bounded. One way to pull the 
system back into physical realm is to couple the vibra- 
tional mode with a thermal bath, see Appendix A. As 
explained above for the TLS-mode case, the following 
damping rate is a measure for the junction stability. 



A negative value indicates on junction instability, as un- 
controlled heating of the mode takes place. Using steady- 
state populations, we proceed and derive the charge cur- 
rent expression, valid only if fc(j < fc^, 
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(67) 



(68) 



Here, Pgs is a vector of steady-state population given 
by Eq. (66). Equations (59)- (68) can be generalized 
to include the interaction of the harmonic mode with 
a dissipative-thermal phonon bath. Appendix A exem- 
plifies this procedure for the anharmonic-mode model. 
In practice, the electronic induced rates fc^ and k'^ in 
Eq. (66) are augmented by a phononic contribution, 
Tph{uJo)[nph{'^o) + 1] and rph(a;o)np/i(wo), respectively. 

Fig. 9 displays the charge current for zero and finite 
Tph strength. In the absence of coupling to the phonon 
bath, the current diverges around A/i = 0.6, where insta- 
bility occurs. For finite Fp/j, the current is larger when 
the vibrational mode is harmonic, compared to the TLS- 
mode case (dotted) , as the electronic energy can be used 
to excite multiple transitions. 

Fig. 10 further demonstrates the "stabilizing" effect 
the interaction with a heat bath has on the harmonic 
mode. We display the population for the states n = 
to n = 3, top to bottom. For Tph = the data is 
presented up to A/x = 0.6, where the population becomes 
unphysical (levels' population goes to zero there since 
infinite number of vibrational states are occupied). This 
point is indicated by an arrow. 

The cooling and heating behavior depicted in Fig. 8 
for the TLS model could be repeated for the present 
harmonic-mode case as well, to yield the same behav- 
ior. The reason is that /3e// is determined by the ra- 
tio of rates, and this ratio is identical in the two mod- 
els. Our conclusion here is thus that including modc- 
anharmonicity is important for transport calculations: 
Harmonic-mode model can lead to unphysical results 
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FIG. 9: Charge current in the harmonic-mode model for 
Tp/j = (dashed line) and Fp^ = 0.05 (full line). For compar- 
ison, we also present the current in the TLS-mode model with 
Fph = 0.05 (dotted line), fiv ~ Pph = 40, other parameters 
are the same as in Fig. 3. 
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region codes instability zones, with negative K^b- These 
diagrams hold for both TLS and harmonic mode cases. 
A reentrant behavior is observed in Fig. 11: For a fixed 
value of Tph, say Tph — 0.005, the junction is stable for 
A/i < 0.6, unstable around 0.6 < Afi < 1.2, while beyond 
that, the junction is operative again. The reason for this 
behavior is that electronic-induced excitation and relax- 
ation rates, fc^ and k^, become both small at large bias, 
thus the thermal bath-induced rates dominate the mode 
dynamics, leading to a normal-thermal like behavior. 

The coupling of the D and A molecular states to the 
metal leads may be further tuned experimentally. In Fig. 
12 we show a stability map of the junction as a function of 
voltage bias and metal-molecule hybridization F^ = F^. 
We explore two situations: (a) The molecular mode is 
perfectly isolated from other vibrations, and (b) Tph is 
finite. In both cases, once metal-molecule coupling is 
large enough, a stable operation sustains. This result 
seems initially counterintuitive, as one expects strongly 
coupled molecules to support high charge and energy cur- 
rents, potentially leading to junction rupture. However, 
the key factor in the formation of vibrational instability 
here is the fact that at certain voltages the vibrational 
excitation rate k^ exceeds the relaxation rate k^. In- 
specting the rates (25), one can analytically prove that 
if the effective density of states is energy independent, 
Ji^{e) = C, which is the case at strong metal-molecule 
coupling, then k^ — fc^ ex C^wq, a positive number. The 
key factor in instability build-up is thus the usage of elec- 
tronic reservoirs with effective DOS [Eq. (24)] peaked 
around different energies, the D and A levels. 



FIG. 10: Population of the first four levels of the harmonic 
mode. Full lines: weak interaction with the heat bath, Tph~ 
10~^. The population becomes unphysical (negative) for 
Afi > 0.6. Dashed lines: strong interaction with heat bath, 
Fph =0.1, lifts the instability. /?,^ — l3ph ~ 40, other parame- 
ters are the same as in Fig. 3. 



(e.g., current divergence) since there is no saturation sit- 
uation for the vibrational mode. In particular, including 
molecular anharmonic aspects is essential for obtaining 
reliable results when simulating junction behavior close 
to the critical bias, where an instability occurs. 
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V. SUMMARY 

We have studied vibrational cooling, heating, and in- 
stability formation in a phonon assisted D-A electron rec- 
tifier junction using a full-counting statistics approach. 
Variants of the basic model were constructed, assuming 
either harmonic or an anharmonic vibrational mode, fur- 
ther allowing energy dissipation to a phononic thermal 
environment. 

Putting together our observations, we present in fig- 
ures 11 and 12 stability maps for the system; the dark 



FIG. 11: Stability diagram. The dark island and the narrow 
strip (at negative bias) are the parametric region in which the 
junction becomes unstable. Other parameters are the same 
as in Fig. 3, besides the temperatures, I3l = Pr = Pph = 40. 

Concluding our observations: (i) We confirmed the 
steady-state entropy production FT for the different 
model variants. This is a non-trivial task since charge 
and energy currents here are not tightly coupled, a result 
of the strong metal-D and A-metal couplings. There- 
fore, one needs to separately count particle number and 
energy transfer in the system, (ii) We derived simple 
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FIG. 12: (a): Stability diagram for Fph = 0. (b); Stability 
diagram with Fph ~ 0.005. The dark region is the paramet- 
ric region in which the junction becomes unstable, K^it < 0. 
Other parameters are the same as in Fig. 3, besides the tem- 
peratures Pl = Pr. = Pph = 40. 



analytical expressions for the charge current, assuming 
either harmonic or an anharmonic vibrational mode. As 
expected, harmonic-mode junctions better conduct since 
the electronic energy can be used to excite multiple vi- 
brational states. An anharmonic mode quickly reaches 
saturation, (iii) We defined an effective temperature for 
the vibration and demonstrated bias induced cooling at 
low bias, A/i < 0.2. (iv) At intermediate voltage bias, 
0.2 < A^ < 0.6, heating effects were observed, and the 
mode effective temperature exceeds the environmental 
temperature, (v) Once the donor and acceptor levels 
switch position for A/z > 0.6, ed > £„, instability de- 
velops: The mode excitation rate exceeds the relaxation 
rate, and the vibrational mode uncontrollably heats, (vi) 
Coupling the vibrational mode to an external thermal 
bath stabilizes the junction. In particular, assuming a 
weak interaction to a phonon bath, Tph — 0.005, junc- 
tion instability is removed for A/i > 1.2; the electronic- 
induced rates diminish and the mode dynamics is con- 
trolled by the thermal bath, (vii) The appearance of 
vibrational instability can be traced down to the metals' 
energy dependent DOS, different at the two ends. 

The simple models described here elucidate the role 
of different factors on vibrational cooling, heating, and 
instability build-up in a D-A electronic rectifier. The 
effects of mode frequency, its interaction with other 
modes, the reservoirs' temperature, metal-molecule cou- 
pling strength, and bias voltage, were examined. While 
the focus of this work has been on vibrational effects, 
the theory developed here could be used for describing 
the coupling of an electronic junction to a cavity mode- 
electromagnetic environment. One could thus reformu- 
late this study and describe cooling, heating and diode- 
like effects in photonic heat conduction [48-52]. Future 



work will be devoted to the study of noise processes in 
phonon-assisted tunneling junctions, with the motiva- 
tion to expose mode properties (harmonicity) through 
the noise characteristics. 
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Appendix A: Full counting statistics for charge and 

energy in the dissipative anharmonic-mode rectifier 

model 



We describe here the derivation of the generating func- 
tion and the charge current for the anharmonic-mode 
bath-coupled rectifier model. 



Ha+b = -z-cTz+o-x (Fe 



J2^ialai + J2 



Fh 



Ct" (Jjj,(Jjf 



J2^c.bib^.{Al) 



Here, the TLS excitation and relaxation processes are 
coupled to both electronic transitions in the junction and 
to a thermal-phononic reservoir, 

Fe = Ky^{X*Xrajar + \*\ialai). 

l,r 

Fb - ^t;a(6l + 6„). (A2) 

Q 

ai , Or are fcrmionic annihilation operators corresponding 
to the left and right metals; the coefficients A^ were de- 
fined in Eq. (8), n denotes the D-A tunneling strength, ha 
is a bosonic operator, describing the a reservoir mode, Va 
quantifies the TLS-bath interaction strength. For more 
details, see text around Eq. (55). 

The impurity (TLS) dynamics can be obtained by us- 
ing a master-equation approach [43]. The procedure in- 
volves a second order perturbation theory in the impurity 
coupling to both electronic and phononic reservoirs. In 
the markovian limit, under the rotating wave approxi- 
mation, one standardly achieves kinetic equations that 
separately account for the electronic (e) and phononic 
(b) relaxation pathways. 



PI - - (fci^o + fcLo) Pi + {K- 
pi+pn = 1- 

The relaxation terms are given by 



IPo 



(A3) 



o«(«»i-«„')t 



'^^Fe{r)FM)dr 



^'^^Fb{T)FbmdT. (A4) 
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Here, e„ is the energy of the nth vibrational level. Elec- 
tron induced rate constants are detailed through Eqs. 
(19)-(25). The thermal bath induced rates can be simi- 
larly put together, 



fci^o = rp;i(a;o)K/i(^o) + 1], 
k'' - k'' 



-'»Jo/3ph 



(A5) 



^ph 



(oj) = [e'^p'i" — 1]^^ denotes the Bose-Einstein dis- 
tribution function and Tph{uj) = ^tt^^^ w^(5(a;a — ui). 
For brevity, we ignore below the direct reference to fre- 
quency. We now define the probability distribution func- 



tion Pt(n, A^, wl, Wfl, gwo), as the probability to find the 
system at time t in state n = 0, 1, with N electrons trans- 
ferred to the right bath, uii, excess energy accumulated 
at the ly bath (i^ = L,R), and qujo energy attained by 
the phonon bath, due to the transfer of q quantas from 
the TLS to this bath. Note that charge conservation be- 
tween the L and R baths is enforced, allowing us to work 
with a single counting field for describing charge transfer 
processes in the steady-state limit. 

We now resolve the associated master equation for the 
two-state population, to its charge and energy contribu- 
tions. The resulting equations are analogous to Eq. (27), 



J 



Pti'^,N,ujL,ujji,quJo) = -Vtil,N,ujL,ujR,qujQ)[kl^Q + TphioJa)[nph{ujQ) + 1]] 

+ / Vt{0,N -l,ujL + e,ujR- e + uJo,quJo)fL{e)['i-- fR{e-ujo)]JLie)JR{e-oJo)de 



+ / Vt{0,N + l,LOL - e + uJo,LOR + e,quJo)fR{e)[l - /^(e - wo)] Jfl(e) ^^(e - wo)de 

J — OO 

Pt{0,N,ujL,ujR,quJo) = -'Pt{0,N,ujL,UJR,quJa)[kQ^^+Tph{uJa)nph{uJo)] 

Pt{i,N - 1,UJL + e,ujR- e-uJa,qujQ)fL{e)[l - fR{e + ujQ)]JL{e)jR{e + ujQ)de 



+ I 7^4(1, A^ + 1,CJL - e- Wo, cj_R + e, (7CJo)/i?(e)[l - /L(e + wo)]J_R(e)-^L(e + wo)& 

-I- Pt{i,N,ujL,ujR,{q-l)ujo)'rph{(^o)[nph{(^o) + M- (A6) 



r 



We Fourier transform this system with respect to charge depends on the energy counting fields i]^ and ^, and the 

charge counting field x- It is a vector with two entries, 
as in Eq. (28), and it satisfies the differential equation 



and energy, 

T^t{n,x,VL,VR,0 



J2 e*^'^ Yl ^"^'^"^ I ^'""^"^ ^'^i 



e 



d\Z) 
dt 



= -W\Z). 



(A8) 



X / e'"«''«dwflPt(n,7V,WL,Wfl,gc^o) (A7) 



to obtain the characteristic function Z{x,TjL,TlR,(,Tt). It The matrix W has the following entries 
I 



W 



"^0— >1 "■" I phlT'ph 



~e'^F^{f^L,m) - e-'^'F+ir^L^VR) - ^ph{nph + l)e*«"" 



k'l^Q +Tphinph + 1) 



(A9) 



where 



and 



Ft{VL,VR.) = / 



^-ieriL ^i(£T'^o)VR 



F^{VL,VR.) = / 

J — ( 



giJJL(e±tJo)g-ier;H 



X /L(e)[l - fR{eTuJo)]JL{e)JRieT(^o)de 

(AlO) 



X [1 - /^(e ± uJo)]fR{e)JL{e ± uJo)jR{e)de. 

(All) 
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The CGF is expressed in terms of the characteristic func- 
tion \Z) a,s 

G(x,'?L,^/?,C)= lim 7ln(/|Z), (A12) 

Practically, it is reached by the negative of the smallest 
eigenvalue of the matrix W, 



G{x,VL,VR,0 



Wll +W2A 



(A13) 



^(Wl,l - W2.2)^ + 4U71,2W2,1 



where the following holds when the counting fields are all 
set to zero, wi^2 — — W2,2, W2.1 = -w^i.i, dwi^2/d{ix)\o = 
fcf-^o^-fcf^o^ a.iiddw2,i/d{ix)\o = k^^i-k^^^. RecaU 
that the rates fcJ^I^J^; are electron-hole generation assisted 
rates, see Eq. (25). We can now identify the levels pop- 
ulation, 



Wis 



-L ph'^ph 



Wii + W2.2 



"■0- 
Pl, 



+ ^1^0 + rph[2nph + 1] 



(A15) 



with Wij the matrix elements in Eq. (A9). The charge 
current is obtained by taking the first derivative with 
respect to (ix). 



(le) 



dG 

dix 



(Wl,l + W2,2) ^ 



W2,l- 



dwi^2 
dix 



wi 2- 



dW2,l 



dix 



(A14) 



X^V,^,i=o 



and similarly for pq = 1 — pi, resulting in the expression 
for the charge current 



(/e>=Pi(A:f::o''-fcf:?o^)+Po(fco' 



L^R 



„fl^i^ 



(A16) 



This result is formally identical to Eq. (48), with the only 
difference that the steady-state TLS population is now 
modified, to include phonon-bath assisted transitions. 
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